1 Setup

1.1 Libs

library(car) 
## Loading required package: carData
library("Rmisc")
## Loading required package: lattice
## Loading required package: plyr
library("ggpubr")
## Loading required package: ggplot2
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
## 
##     mutate
library("cowplot")
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
# library("lme4")
library("phyloseq")
library("glmmTMB")
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.9.11
## Current TMB version is 1.9.14
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
library("emmeans")
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
# library(multcomp)
# #install.packages("multcompView")
# library(multcompView)
library("bbmle")
## Loading required package: stats4
library("DHARMa")
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
# library("nlme")
#install.packages("bestNormalize")
library("bestNormalize")

1.2 Data

setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/mosquito_microbes/04.microbiome_analysis/04a.diversity")

ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.clean 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 161 taxa and 458 samples ]
## sample_data() Sample Data:       [ 458 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 161 taxa by 11 taxonomic ranks ]
samdf <- data.frame(ps.clean@sam_data)

##made below
#df.div <- read.csv("div.trim.csv",row.names=1)

2 Alpha diversity calculations

df.all <- data.frame(estimate_richness(ps.clean, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))

df.all$sample_name <- rownames(df.all)

df.all.div <- merge(df.all,samdf,by="sample_name") #add sample data

#shannon diversity divided by species richness
df.all.div$even <- df.all.div$Shannon/(log(df.all.div$Observed))

df.div <- df.all.div

#adding most recent sample counts
sums <- data.frame(sample_sums(ps.clean))
colnames(sums) <- c("lib_size_trim")
sums$sample_name <- row.names(sums)

df.div <- merge(df.div,sums,by="sample_name")

#write.csv(df.div,"div.trim.csv")

3 Analyzing alpha metrics

3.1 All data

Very strange/interesting pattern where the water types that shouldn’t have a lot of diversity have a lot of OTUs, but not a lot of reads…

ggplot(df.div,aes(x=type,y=Observed,color=infusion))+
  geom_boxplot()+
  geom_jitter()

ggplot(df.div,aes(x=type,y=lib_size_trim,color=infusion))+
  geom_boxplot()+
  geom_jitter(position=position_jitterdodge())

3.2 Plotz

#just experimental types
df.div.exp <- subset(df.div,type=="Microbial Water"|type=="A.albopictus")

df.div.exp$inf.temp <- paste0(df.div.exp$infusion,df.div.exp$temperature)

df.div.exp$type <- sub("A.albopictus","Mosquitoes",df.div.exp$type)
df.div.exp$type <- sub("Microbial Water","Meso. water",df.div.exp$type)

df.div.exp$type <- factor(df.div.exp$type,levels=c("Mosquitoes","Meso. water"))

3.2.1 Richness

df.div.exp.se <- summarySE(df.div.exp,measurevar="Observed",groupvars=c("inf.temp","infusion","temperature","type"))

##dots & error bars
gg.rich <- ggplot(df.div.exp,aes(x=infusion,y=Observed,fill=inf.temp,shape=inf.temp))+
  geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  ylim(4,34)+
  facet_wrap(~type)+
  theme_cowplot()+
  ylab("ASV richness")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))
gg.rich

##water by time
df.div.water <- subset(df.div.exp,type=="Meso. water")
gg.mw.rich.time <- ggplot(df.div.water,aes(x=infusion,y=Observed,fill=inf.temp))+
  geom_boxplot()+
  #geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  #geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  #geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~day)+
  theme_cowplot()+
  ylab("ASV richness")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  ggtitle("Mesocosm water")
gg.mw.rich.time

3.2.2 Simpson’s

df.div.exp.se <- summarySE(df.div.exp,measurevar="InvSimpson",groupvars=c("inf.temp","infusion","temperature","type"))

##dots & error bars
gg.simp <- ggplot(df.div.exp,aes(x=infusion,y=InvSimpson,fill=inf.temp,shape=inf.temp))+
  geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  geom_errorbar(data=df.div.exp.se,aes(ymax=InvSimpson+se,ymin=InvSimpson-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~type)+
  theme_cowplot()+
  ylab("Simpson's (inv.)")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))+
  scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL, cool","OL, warm","SG, cool","SG, warm","PW, cool","PW, warm"))
gg.simp

gg.mw.simp.time <- ggplot(df.div.water,aes(x=infusion,y=InvSimpson,fill=inf.temp))+
  geom_boxplot()+
  #geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
  #geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
  #geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
  #geom_boxplot(outlier.shape=NA,alpha=0.5)+
  facet_wrap(~day)+
  theme_cowplot()+
  ylab("Simpson's index (inv.)")+
  xlab("Infusion")+
  scale_x_discrete(labels=c("OL","SG","PW"))+
  scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
  ggtitle("Mesocosm water")

3.2.3 Multi-panel (Fig 2 subpanel & Fig S10)

gg.div <- ggarrange(gg.rich,gg.simp,common.legend=T,legend="right",labels=c("a","b"))
gg.div

##just water things
ggarrange(gg.mw.rich.time,gg.mw.simp.time,common.legend=T,legend="right",nrow=2)

#ggsave("div.water.time.pdf",width=6,height=6)

3.3 Statz

3.3.1 Richness

3.3.1.1 Water richness stats

df.div.mw <- subset(df.div.exp,type=="Meso. water")

df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)
## 'data.frame':    211 obs. of  30 variables:
##  $ sample_name      : chr  "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
##  $ Observed         : num  25 28 22 28 31 31 28 29 27 21 ...
##  $ Shannon          : num  1.63 1.53 1.58 1.53 1.75 ...
##  $ InvSimpson       : num  2.67 2.4 2.44 2.32 2.93 ...
##  $ orgname          : chr  "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
##  $ mesocosm         : chr  "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
##  $ type             : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
##  $ stage            : chr  "Water" "Water" "Water" "Water" ...
##  $ sex              : chr  "" "" "" "" ...
##  $ date.collected   : chr  "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
##  $ day              : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
##  $ abdomen.length.mm: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ wing.length      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ hindleg.length   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ special          : chr  "" "" "" "" ...
##  $ dispersal        : chr  "N" "N" "D" "D" ...
##  $ temperature      : chr  "C" "C" "C" "C" ...
##  $ infusion         : chr  "OL" "OL" "OL" "OL" ...
##  $ raw_miseq_reads  : int  30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
##  $ X16scq           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct1         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Act_ct           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ notes            : chr  "" "" "" "" ...
##  $ newday           : chr  "early" "early" "early" "early" ...
##  $ lib_size         : num  25797 34139 15861 49298 44559 ...
##  $ is.neg           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ even             : num  0.506 0.458 0.511 0.46 0.51 ...
##  $ lib_size_trim    : num  25797 34134 15861 49290 44524 ...
##  $ inf.temp         : chr  "OLC" "OLC" "OLC" "OLC" ...
#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
#shapiro.test(df.div.mw$Observed)

wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)

AICtab(wrich1,wrich2,wrich3)
##        dAIC  df
## wrich3   0.0 9 
## wrich2 302.7 9 
## wrich1 312.9 8
##Gaussian best

##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus temperature 
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus time 
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)

anova(wrich1,wrich1nd) #0.3931
## Data: df.div.mw
## Models:
## wrich1nd: Observed ~ day + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nd:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1nd  8 1035.4 1062.2 -509.69   1019.4                         
## wrich1    9 1036.7 1066.8 -509.33   1018.7 0.7293      1     0.3931
anova(wrich1,wrich1ni) #2.2e-16***
## Data: df.div.mw
## Models:
## wrich1ni: Observed ~ day + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1ni:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1ni  7 1118.0 1141.5 -552.02   1104.0                             
## wrich1    9 1036.7 1066.8 -509.33   1018.7 85.373      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrich1,wrich1nt) #0.4119
## Data: df.div.mw
## Models:
## wrich1nt: Observed ~ day + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nt:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1nt  8 1035.3 1062.2 -509.67   1019.3                         
## wrich1    9 1036.7 1066.8 -509.33   1018.7 0.6732      1     0.4119
anova(wrich1,wrich1nnd) #0.001874**
## Data: df.div.mw
## Models:
## wrich1nnd: Observed ~ temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nnd:     (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## wrich1nnd  7 1045.2 1068.7 -515.61   1031.2                            
## wrich1     9 1036.7 1066.8 -509.33   1018.7 12.559      2   0.001874 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day          12.9404  2   0.001549 ** 
## temperature   0.6742  1   0.411578    
## infusion    151.4172  2  < 2.2e-16 ***
## dispersal     0.7306  1   0.392704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: Observed
#                Chisq Df Pr(>Chisq)    
# day          12.9404  2   0.001549 ** 
# temperature   0.6742  1   0.411577    
# infusion    151.4172  2  < 2.2e-16 ***
# dispersal     0.7306  1   0.392704    
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full

##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)

anova(wrich1,wrich1.int) #2.86e-10***
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int:     (1 | mesocosm), zi=~0, disp=~1
##            Df     AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1      9 1036.66 1066.8 -509.33  1018.66                             
## wrich1.int 13  994.79 1038.4 -484.40   968.79 49.863      4  3.857e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrich1,wrich1.int1) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int1: Observed ~ day + infusion * temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int1:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int1 11 1037.2 1074.1 -507.62   1015.2 3.4154      2     0.1813
anova(wrich1,wrich1.int2) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int2: Observed ~ day + infusion + temperature * dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int2:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int2 10 1037.1 1070.7 -508.57   1017.1 1.5253      1     0.2168
anova(wrich1,wrich1.int3) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int3: Observed ~ day * temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int3:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int3 11 1039.0 1075.9 -508.52   1017.0 1.6226      2     0.4443
anova(wrich1,wrich1.int4) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int4: Observed ~ day * dispersal + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int4:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int4 11 1039.3 1076.2 -508.66   1017.3 1.3341      2     0.5132
anova(wrich1,wrich1.int5) #ns
## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int5: Observed ~ day + dispersal * infusion + temperature + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int5:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1036.7 1066.8 -509.33   1018.7                         
## wrich1.int5 11 1037.8 1074.7 -507.90   1015.8 2.8599      2     0.2393
Anova(wrich1.int)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                 Chisq Df Pr(>Chisq)    
## day           17.4815  2  0.0001599 ***
## infusion     167.8498  2  < 2.2e-16 ***
## temperature    0.7595  1  0.3834697    
## dispersal      0.7626  1  0.3825074    
## day:infusion  59.0003  4  4.705e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                 Chisq Df Pr(>Chisq)    
# day           17.4811  2    0.00016 ***
# infusion     167.8537  2  < 2.2e-16 ***
# temperature    0.7595  1    0.38348    
# dispersal      0.7627  1    0.38248    
# day:infusion  58.9948  4  4.717e-12 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)

wrich1.int.notem <- glmmTMB(Observed~day*infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1.int,wrich1.int.notem) #0.38
## Data: df.div.mw
## Models:
## wrich1.int.notem: Observed ~ day * infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int.notem:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int:     (1 | mesocosm), zi=~0, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1.int.notem 12 993.55 1033.8 -484.78   969.55                         
## wrich1.int       13 994.79 1038.4 -484.40   968.79 0.7567      1     0.3844
wrich1.int.nodis <- glmmTMB(Observed~day*infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1.int,wrich1.int.nodis) #0.38
## Data: df.div.mw
## Models:
## wrich1.int.nodis: Observed ~ day * infusion + temperature + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int.nodis:     (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int:     (1 | mesocosm), zi=~0, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1.int.nodis 12 993.55 1033.8 -484.78   969.55                         
## wrich1.int       13 994.79 1038.4 -484.40   968.79 0.7574      1     0.3841
##model checking
shapiro.test(residuals(wrich1.int)) #awesome
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wrich1.int)
## W = 0.99201, p-value = 0.3043
qqnorm(residuals(wrich1.int))
qqline(residuals(wrich1.int), col="red")

plotResiduals(wrich1.int)

wrich1.int.resid <- simulateResiduals(fittedModel = wrich1.int, plot = T)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)

multcomp::cld(wrich1.int.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  SW       C             19.1 0.351 198     18.5     19.8  1    
##  SW       H             19.5 0.353 198     18.8     20.1  1    
##  SG       C             21.2 0.355 198     20.5     21.9   2   
##  SG       H             21.5 0.357 198     20.8     22.2   2   
##  OL       C             24.7 0.351 198     24.0     25.4    3  
##  OL       H             25.0 0.353 198     24.3     25.7    3  
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean    SE  df lower.CL upper.CL .group
 # SW       C             19.1 0.351 198     18.5     19.8  1    
 # SW       H             19.5 0.353 198     18.8     20.1  1    
 # SG       C             21.2 0.355 198     20.5     21.9   2   
 # SG       H             21.5 0.357 198     20.8     22.2   2   
 # OL       C             24.7 0.351 198     24.0     25.4    3  
 # OL       H             25.0 0.353 198     24.3     25.7    3  

3.3.1.2 Water richness stats - rarefied

ps.rare <- readRDS("../../02.process_asvs/ps.clean.trim.rare9200.rds")
ps.rare
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 102 taxa and 406 samples ]
## sample_data() Sample Data:       [ 406 samples by 24 sample variables ]
## tax_table()   Taxonomy Table:    [ 102 taxa by 11 taxonomic ranks ]
samdf.rare <- data.frame(ps.rare@sam_data)

df.rare <- data.frame(estimate_richness(ps.rare, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))

df.rare$sample_name <- rownames(df.rare)

df.rare.div <- merge(df.rare,samdf.rare,by="sample_name") #add sample data

df.rare.div.exp <- subset(df.rare.div,type=="Microbial Water"|type=="A.albopictus")

Same as above but without offset for library size

df.rare.div.mw <- subset(df.rare.div.exp,type=="Microbial Water")

df.rare.div.mw$day <- as.factor(df.rare.div.mw$day)
str(df.rare.div.mw)
## 'data.frame':    211 obs. of  27 variables:
##  $ sample_name      : chr  "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
##  $ Observed         : num  25 27 22 28 31 31 28 29 27 21 ...
##  $ Shannon          : num  1.64 1.51 1.56 1.5 1.77 ...
##  $ InvSimpson       : num  2.69 2.36 2.41 2.27 2.98 ...
##  $ orgname          : chr  "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
##  $ mesocosm         : chr  "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
##  $ type             : chr  "Microbial Water" "Microbial Water" "Microbial Water" "Microbial Water" ...
##  $ stage            : chr  "Water" "Water" "Water" "Water" ...
##  $ sex              : chr  "" "" "" "" ...
##  $ date.collected   : chr  "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
##  $ day              : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
##  $ abdomen.length.mm: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ wing.length      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ hindleg.length   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ special          : chr  "" "" "" "" ...
##  $ dispersal        : chr  "N" "N" "D" "D" ...
##  $ temperature      : chr  "C" "C" "C" "C" ...
##  $ infusion         : chr  "OL" "OL" "OL" "OL" ...
##  $ raw_miseq_reads  : int  30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
##  $ X16scq           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct1         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Act_ct           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ notes            : chr  "" "" "" "" ...
##  $ newday           : chr  "early" "early" "early" "early" ...
##  $ lib_size         : num  25797 34139 15861 49298 44559 ...
##  $ is.neg           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
#hist(df.rare.div.mw$Observed)
#shapiro.test(log(df.rare.div.mw$Observed))
#shapiro.test(df.rare.div.mw$Observed)

#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="poisson",data=df.rare.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.rare.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)

#AICtab(wrich1,wrich2,wrich3)

##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+(1|mesocosm), data=df.rare.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus temperature 
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus time 
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)

anova(wrich1,wrich1nd) #0.466
## Data: df.rare.div.mw
## Models:
## wrich1nd: Observed ~ day + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1nd  8 1066.8 1093.7 -525.43   1050.8                         
## wrich1    9 1068.3 1098.5 -525.16   1050.3 0.5315      1      0.466
anova(wrich1,wrich1ni) #2.2e-16***
## Data: df.rare.div.mw
## Models:
## wrich1ni: Observed ~ day + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1ni  7 1145.2 1168.7 -565.62   1131.2                             
## wrich1    9 1068.3 1098.5 -525.16   1050.3 80.917      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrich1,wrich1nt) #0.4789
## Data: df.rare.div.mw
## Models:
## wrich1nt: Observed ~ day + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1nt  8 1066.8 1093.6 -525.41   1050.8                         
## wrich1    9 1068.3 1098.5 -525.16   1050.3 0.5014      1     0.4789
anova(wrich1,wrich1nnd) #0.006693**
## Data: df.rare.div.mw
## Models:
## wrich1nnd: Observed ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## wrich1nnd  7 1074.3 1097.8 -530.17   1060.3                            
## wrich1     9 1068.3 1098.5 -525.16   1050.3 10.013      2   0.006693 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                Chisq Df Pr(>Chisq)    
## day          10.2547  2   0.005932 ** 
## temperature   0.5020  1   0.478603    
## infusion    133.7700  2  < 2.2e-16 ***
## dispersal     0.5322  1   0.465696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                Chisq Df Pr(>Chisq)    
# day          10.2547  2   0.005932 ** 
# temperature   0.5020  1   0.478604    
# infusion    133.7700  2  < 2.2e-16 ***
# dispersal     0.5322  1   0.465695    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full

##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+(1|mesocosm),data=df.rare.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+(1|mesocosm),data=df.rare.div.mw)

anova(wrich1,wrich1.int) #4.35e-11***
## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wrich1      9 1068.3 1098.5 -525.16  1050.32                             
## wrich1.int 13 1021.9 1065.5 -497.96   995.93 54.391      4  4.358e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wrich1,wrich1.int1) #ns
## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int1: Observed ~ day + infusion * temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1068.3 1098.5 -525.16   1050.3                         
## wrich1.int1 11 1068.7 1105.5 -523.33   1046.7 3.6606      2     0.1604
anova(wrich1,wrich1.int2) #ns
## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int2: Observed ~ day + infusion + temperature * dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1       9 1068.3 1098.5 -525.16   1050.3                        
## wrich1.int2 10 1068.5 1102.0 -524.24   1048.5 1.836      1     0.1754
anova(wrich1,wrich1.int3) #ns
## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int3: Observed ~ day * temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1068.3 1098.5 -525.16   1050.3                         
## wrich1.int3 11 1071.0 1107.8 -524.48   1049.0 1.3561      2     0.5076
anova(wrich1,wrich1.int4) #ns
## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int4: Observed ~ day * dispersal + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1068.3 1098.5 -525.16   1050.3                         
## wrich1.int4 11 1071.5 1108.3 -524.74   1049.5 0.8411      2     0.6567
anova(wrich1,wrich1.int5) #ns
## Data: df.rare.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int5: Observed ~ day + dispersal * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1       9 1068.3 1098.5 -525.16   1050.3                         
## wrich1.int5 11 1070.2 1107.1 -524.10   1048.2 2.1147      2     0.3474
Anova(wrich1.int)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                 Chisq Df Pr(>Chisq)    
## day           13.9846  2  0.0009189 ***
## infusion     154.5864  2  < 2.2e-16 ***
## temperature    0.5795  1  0.4464952    
## dispersal      0.5809  1  0.4459485    
## day:infusion  64.8561  4  2.759e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#                 Chisq Df Pr(>Chisq)    
# day           13.9847  2  0.0009189 ***
# infusion     154.5860  2  < 2.2e-16 ***
# temperature    0.5795  1  0.4464988    
# dispersal      0.5809  1  0.4459484    
# day:infusion  64.8562  4  2.759e-13 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

wrich1.int.notem <- glmmTMB(Observed~day*infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
anova(wrich1.int,wrich1.int.notem) #0.45
## Data: df.rare.div.mw
## Models:
## wrich1.int.notem: Observed ~ day * infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1.int.notem 12 1020.5 1060.7 -498.25   996.51                         
## wrich1.int       13 1021.9 1065.5 -497.96   995.93 0.5777      1     0.4472
wrich1.int.nodis <- glmmTMB(Observed~day*infusion+temperature+(1|mesocosm),data=df.rare.div.mw)
anova(wrich1.int,wrich1.int.nodis) #0.45
## Data: df.rare.div.mw
## Models:
## wrich1.int.nodis: Observed ~ day * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wrich1.int.nodis 12 1020.5 1060.7 -498.25   996.51                         
## wrich1.int       13 1021.9 1065.5 -497.96   995.93 0.5778      1     0.4472
#tapply(df.rare.div.mw$log(Observed), df.rare.div.mw$infusion, mean)
#tapply(df.rare.div.mw$log(Observed), df.rare.div.mw$day, mean)

##model checking
shapiro.test(residuals(wrich1.int)) #awesome
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wrich1.int)
## W = 0.99017, p-value = 0.1609
qqnorm(residuals(wrich1.int))
qqline(residuals(wrich1.int), col="red")

plotResiduals(wrich1.int)

wrich1.int.resid <- simulateResiduals(fittedModel = wrich1.int, plot = T)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)

multcomp::cld(wrich1.int.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  SW       C             19.1 0.371 198     18.3     19.8  1    
##  SW       H             19.4 0.373 198     18.6     20.1  12   
##  SG       C             20.9 0.374 198     20.1     21.6   23  
##  SG       H             21.1 0.377 198     20.4     21.9    3  
##  OL       C             24.6 0.371 198     23.9     25.4     4 
##  OL       H             24.9 0.373 198     24.2     25.7     4 
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

3.3.1.3 Mosquito richness stats

df.div.mq <- subset(df.div.exp,type=="Mosquitoes")

hist(df.div.mq$Observed)

shapiro.test(df.div.mq$Observed)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mq$Observed
## W = 0.89464, p-value = 1.697e-10
##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mq)

AICtab(mrich1,mrich2,mrich3) #gaussian again
##        dAIC df
## mrich1  0.0 10
## mrich2 47.1 10
## mrich3 48.3 9
##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##no dispersal 
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus infusion 
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus temperature 
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus time 
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)

anova(mrich1,mrich1nd) #0.24
## Data: df.div.mq
## Models:
## mrich1nd: Observed ~ newday + temperature + infusion + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1nd:     (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1:     offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1nd  9 917.87 947.32 -449.93   899.87                         
## mrich1   10 918.46 951.19 -449.23   898.46 1.4056      1     0.2358
anova(mrich1,mrich1ni) #0.78
## Data: df.div.mq
## Models:
## mrich1ni: Observed ~ newday + temperature + dispersal + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1ni:     (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1:     offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1ni  8 914.95 941.14 -449.48   898.95                         
## mrich1   10 918.46 951.19 -449.23   898.46 0.4942      2     0.7811
anova(mrich1,mrich1nt) #0.64
## Data: df.div.mq
## Models:
## mrich1nt: Observed ~ newday + infusion + dispersal + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1nt:     (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1:     offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1nt  9 916.68 946.14 -449.34   898.68                         
## mrich1   10 918.46 951.19 -449.23   898.46 0.2243      1     0.6358
anova(mrich1,mrich1nnd) #0.40
## Data: df.div.mq
## Models:
## mrich1nnd: Observed ~ temperature + infusion + dispersal + sex + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1nnd:     (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1:     offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1nnd  8 916.27 942.46 -450.14   900.27                        
## mrich1    10 918.46 951.19 -449.23   898.46 1.814      2     0.4037
anova(mrich1,mrich1ns) #0.53
## Data: df.div.mq
## Models:
## mrich1ns: Observed ~ newday + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## mrich1ns:     (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## mrich1:     offset(log(lib_size_trim)) + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1ns  9 916.86 946.32 -449.43   898.86                         
## mrich1   10 918.46 951.19 -449.23   898.46 0.3989      1     0.5277
Anova(mrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##              Chisq Df Pr(>Chisq)
## newday      1.8273  2     0.4011
## temperature 0.2244  1     0.6357
## infusion    0.4926  2     0.7817
## dispersal   1.4646  1     0.2262
## sex         0.3992  1     0.5275
shapiro.test(residuals(mrich1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mrich1)
## W = 0.90866, p-value = 1.328e-09
qqnorm(residuals(mrich1))
qqline(residuals(mrich1), col="red")

mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T) 

##posthoc things
mrich1.em <- emmeans(mrich1,~infusion*temperature)
multcomp::cld(mrich1.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  OL       H             8.34 0.417 185     7.51     9.16  1    
##  SG       H             8.49 0.447 185     7.61     9.37  1    
##  OL       C             8.56 0.340 185     7.89     9.23  1    
##  SG       C             8.72 0.490 185     7.75     9.68  1    
##  SW       H             8.77 0.479 185     7.82     9.71  1    
##  SW       C             8.99 0.564 185     7.88    10.10  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
##no differences

3.3.1.4 Mosquito richness stats - rarefied

Same as above but without offset for library size

df.rare.div.mq <- subset(df.rare.div.exp,type=="A.albopictus")

#hist(df.rare.div.mq$Observed)
#shapiro.test(log(df.rare.div.mq$Observed))
#shapiro.test(df.rare.div.mq$Observed)

#mrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="poisson",data=df.rare.div.mq)
#mrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.rare.div.mq)
#mrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mq)

#AICtab(mrich1,mrich2,mrich3)

##full model
mrich1<-glmmTMB(Observed~newday+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mq)
##no dispersal
mrich1nd<-glmmTMB(Observed~newday+temperature+infusion+(1|mesocosm), data=df.rare.div.mq)
##minus infusion
mrich1ni<-glmmTMB(Observed~newday+temperature+dispersal+(1|mesocosm), data=df.rare.div.mq)
##minus temperature 
mrich1nt<-glmmTMB(Observed~newday+infusion+dispersal+(1|mesocosm), data=df.rare.div.mq)
##minus time 
mrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mq)

anova(mrich1,mrich1nd) #0.18
## Data: df.rare.div.mq
## Models:
## mrich1nd: Observed ~ newday + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1nd  8 914.40 940.58 -449.2   898.40                         
## mrich1    9 914.61 944.06 -448.3   896.61 1.7942      1     0.1804
anova(mrich1,mrich1ni) #0.93
## Data: df.rare.div.mq
## Models:
## mrich1ni: Observed ~ newday + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## mrich1ni  7 910.75 933.66 -448.37   896.75                        
## mrich1    9 914.61 944.06 -448.30   896.61 0.144      2     0.9305
anova(mrich1,mrich1nt) #0.81
## Data: df.rare.div.mq
## Models:
## mrich1nt: Observed ~ newday + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1nt  8 912.67 938.85 -448.33   896.67                         
## mrich1    9 914.61 944.06 -448.30   896.61 0.0605      1     0.8056
anova(mrich1,mrich1nnd) #0.27
## Data: df.rare.div.mq
## Models:
## mrich1nnd: Observed ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1nnd  7 913.23 936.14 -449.61   899.23                         
## mrich1     9 914.61 944.06 -448.30   896.61 2.6238      2     0.2693
Anova(mrich1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##              Chisq Df Pr(>Chisq)
## newday      2.7059  2     0.2585
## temperature 0.0606  1     0.8056
## infusion    0.1441  2     0.9305
## dispersal   1.8600  1     0.1726
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: Observed
#              Chisq Df Pr(>Chisq)
# newday      2.7059  2     0.2585
# temperature 0.0606  1     0.8056
# infusion    0.1441  2     0.9305
# dispersal   1.8600  1     0.1726

#AICtab(mrich1,mrich1nd,mrich1ni,mrich1nt,mrich1nnd) #no temp best... almost tied with no dispersal, then full

##interactions?
mrich1.int <- glmmTMB(Observed~newday*infusion+temperature+dispersal+(1|mesocosm),data=df.rare.div.mq)
## dropping columns from rank-deficient conditional model: newdaymid:infusionSW
mrich1.int1 <- glmmTMB(Observed~newday+infusion*temperature+dispersal+(1|mesocosm),data=df.rare.div.mq)
mrich1.int2 <- glmmTMB(Observed~newday+infusion+temperature*dispersal+(1|mesocosm),data=df.rare.div.mq)
mrich1.int3 <- glmmTMB(Observed~newday*temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mq)
## dropping columns from rank-deficient conditional model: newdaymid:temperatureH
mrich1.int4 <- glmmTMB(Observed~newday*dispersal+temperature+infusion+(1|mesocosm),data=df.rare.div.mq)
mrich1.int5 <- glmmTMB(Observed~newday+dispersal*infusion+temperature+(1|mesocosm),data=df.rare.div.mq)

anova(mrich1,mrich1.int) #ns
## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
## mrich1.int: Observed ~ newday * infusion + temperature + dispersal + (1 | , zi=~0, disp=~1
## mrich1.int:     mesocosm), zi=~0, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1      9 914.61 944.06 -448.30   896.61                         
## mrich1.int 12 918.16 957.43 -447.08   894.16 2.4499      3     0.4844
anova(mrich1,mrich1.int1) #ns
## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
## mrich1.int1: Observed ~ newday + infusion * temperature + dispersal + (1 | , zi=~0, disp=~1
## mrich1.int1:     mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1       9 914.61 944.06 -448.30   896.61                         
## mrich1.int1 11 918.12 954.12 -448.06   896.12 0.4879      2     0.7835
anova(mrich1,mrich1.int2) #ns
## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
## mrich1.int2: Observed ~ newday + infusion + temperature * dispersal + (1 | , zi=~0, disp=~1
## mrich1.int2:     mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1       9 914.61 944.06 -448.30   896.61                         
## mrich1.int2 10 915.93 948.66 -447.97   895.93 0.6735      1     0.4118
anova(mrich1,mrich1.int3) #ns
## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
## mrich1.int3: Observed ~ newday * temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1.int3:     mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1       9 914.61 944.06 -448.30   896.61                         
## mrich1.int3 10 916.54 949.27 -448.27   896.54 0.0623      1     0.8029
anova(mrich1,mrich1.int4) #ns
## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
## mrich1.int4: Observed ~ newday * dispersal + temperature + infusion + (1 | , zi=~0, disp=~1
## mrich1.int4:     mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1       9 914.61 944.06 -448.30   896.61                         
## mrich1.int4 11 917.74 953.74 -447.87   895.74 0.8642      2     0.6491
anova(mrich1,mrich1.int5) #ns
## Data: df.rare.div.mq
## Models:
## mrich1: Observed ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## mrich1:     mesocosm), zi=~0, disp=~1
## mrich1.int5: Observed ~ newday + dispersal * infusion + temperature + (1 | , zi=~0, disp=~1
## mrich1.int5:     mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mrich1       9 914.61 944.06 -448.30   896.61                         
## mrich1.int5 11 916.84 952.84 -447.42   894.84 1.7644      2     0.4139
Anova(mrich1.int)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Observed
##                  Chisq Df Pr(>Chisq)
## newday          2.7402  2     0.2541
## infusion        0.1459  2     0.9297
## temperature     0.0121  1     0.9125
## dispersal       2.4306  1     0.1190
## newday:infusion 2.4653  3     0.4816
# Response: Observed
#                  Chisq Df Pr(>Chisq)
# newday          2.7401  2     0.2541
# infusion        0.1459  2     0.9296
# temperature     0.0121  1     0.9125
# dispersal       2.4306  1     0.1190
# newday:infusion 2.4653  3     0.4816

#tapply(df.rare.div.mq$log(Observed), df.rare.div.mq$infusion, mean)
#tapply(df.rare.div.mq$log(Observed), df.rare.div.mq$day, mean)

##model checking
shapiro.test(residuals(mrich1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mrich1)
## W = 0.90715, p-value = 1.054e-09
qqnorm(residuals(mrich1))
qqline(residuals(mrich1.int), col="red")

plotResiduals(mrich1)

mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
mrich1.em <- emmeans(mrich1,~infusion+temperature)

multcomp::cld(mrich1.em)
##  infusion temperature emmean    SE  df lower.CL upper.CL .group
##  OL       H             7.91 0.410 186     7.10     8.72  1    
##  SG       H             8.02 0.435 186     7.16     8.88  1    
##  OL       C             8.03 0.328 186     7.38     8.67  1    
##  SW       H             8.13 0.463 186     7.22     9.04  1    
##  SG       C             8.14 0.478 186     7.19     9.08  1    
##  SW       C             8.24 0.552 186     7.16     9.33  1    
## 
## Results are averaged over the levels of: newday, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

3.3.2 Simpson’s

3.3.2.1 Water Simpson’s stats

hist(df.div.mw$InvSimpson)

shapiro.test(df.div.mw$InvSimpson)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$InvSimpson
## W = 0.96756, p-value = 9.003e-05
hist(log(df.div.mw$InvSimpson))

shapiro.test(log(df.div.mw$InvSimpson))
## 
##  Shapiro-Wilk normality test
## 
## data:  log(df.div.mw$InvSimpson)
## W = 0.95673, p-value = 5.166e-06
##normalizing because residuals bad below
bestNormalize(df.div.mw$InvSimpson) #orderNorm
## Best Normalizing transformation with 211 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.5509
##  - Box-Cox: 1.4324
##  - Center+scale: 1.4561
##  - Double Reversed Log_b(x+a): 1.6847
##  - Exp(x): 3.3273
##  - Log_b(x+a): 1.6108
##  - orderNorm (ORQ): 1.1559
##  - sqrt(x + a): 1.4286
##  - Yeo-Johnson: 1.4424
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## 1.359 1.947 2.669 3.225 5.130
df.div.mw$simp.norm <- bestNormalize(df.div.mw$InvSimpson)$x.t
shapiro.test(df.div.mw$simp.norm)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.div.mw$simp.norm
## W = 0.99978, p-value = 1
##just replaced InvSimpson with simp.norm
##full model
wsimp1 <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wsimp1nd <- glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm),data=df.div.mw)
##minus infusion
wsimp1ni <- glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm),data=df.div.mw)
##minus temperature 
wsimp1nt <- glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##minus time 
wsimp1nnd <- glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)

anova(wsimp1,wsimp1nd) #0.86
## Data: df.div.mw
## Models:
## wsimp1nd: simp.norm ~ day + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1nd  8 388.09 414.90 -186.04   372.09                         
## wsimp1    9 390.08 420.24 -186.04   372.08 0.0091      1      0.924
anova(wsimp1,wsimp1ni) #2.2e-16***
## Data: df.div.mw
## Models:
## wsimp1ni: simp.norm ~ day + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1ni  7 489.84 513.30 -237.92   475.84                             
## wsimp1    9 390.08 420.24 -186.04   372.08 103.76      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wsimp1,wsimp1nt) #9.95e-10***
## Data: df.div.mw
## Models:
## wsimp1nt: simp.norm ~ day + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1nt  8 420.66 447.48 -202.33   404.66                             
## wsimp1    9 390.08 420.24 -186.04   372.08 32.585      1  1.141e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wsimp1,wsimp1nnd) #2.65e-09***
## Data: df.div.mw
## Models:
## wsimp1nnd: simp.norm ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1nnd  7 428.81 452.27 -207.40   414.81                             
## wsimp1     9 390.08 420.24 -186.04   372.08 42.732      2  5.259e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wsimp1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                Chisq Df Pr(>Chisq)    
## day          49.7901  2  1.542e-11 ***
## temperature  41.3001  1  1.306e-10 ***
## infusion    229.1056  2  < 2.2e-16 ***
## dispersal     0.0091  1      0.924    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#                Chisq Df Pr(>Chisq)    
# day          45.4454  2  1.354e-10 ***
# temperature  48.9643  1  2.607e-12 ***
# infusion    276.6117  2  < 2.2e-16 ***
# dispersal     0.0291  1     0.8644    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##interaction?
wsimp1.int <- glmmTMB(simp.norm~day*infusion+temperature+dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int1 <- glmmTMB(simp.norm~day+infusion*temperature+dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day+infusion+temperature*dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int3 <- glmmTMB(simp.norm~day*temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
wsimp1.int4 <- glmmTMB(simp.norm~day*dispersal+temperature+infusion+(1|mesocosm),data=df.div.mw)
wsimp1.int5 <- glmmTMB(simp.norm~day+dispersal*infusion+temperature+(1|mesocosm),data=df.div.mw)

anova(wsimp1,wsimp1.int) #4.27e-09***
## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1      9 390.08 420.24 -186.04   372.08                             
## wsimp1.int 13 349.46 393.04 -161.73   323.46 48.615      4  7.026e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wsimp1,wsimp1.int1) #0.026* #goes away with rarefied data
## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int1: simp.norm ~ day + infusion * temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 390.08 420.24 -186.04   372.08                         
## wsimp1.int1 11 389.97 426.84 -183.98   367.97 4.1096      2     0.1281
anova(wsimp1,wsimp1.int2) #ns
## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int2: simp.norm ~ day + infusion + temperature * dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 390.08 420.24 -186.04   372.08                         
## wsimp1.int2 10 391.39 424.90 -185.69   371.39 0.6898      1     0.4062
anova(wsimp1,wsimp1.int3) #ns
## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int3: simp.norm ~ day * temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 390.08 420.24 -186.04   372.08                         
## wsimp1.int3 11 393.77 430.64 -185.89   371.77 0.3052      2     0.8585
anova(wsimp1,wsimp1.int4) #ns
## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int4: simp.norm ~ day * dispersal + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 390.08 420.24 -186.04   372.08                         
## wsimp1.int4 11 392.09 428.96 -185.04   370.09 1.9908      2     0.3696
anova(wsimp1,wsimp1.int5) #ns
## Data: df.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int5: simp.norm ~ day + dispersal * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 390.08 420.24 -186.04   372.08                         
## wsimp1.int5 11 393.48 430.35 -185.74   371.48 0.5976      2     0.7417
Anova(wsimp1.int)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                 Chisq Df Pr(>Chisq)    
## day           70.7469  2  4.340e-16 ***
## infusion     227.4983  2  < 2.2e-16 ***
## temperature   39.6143  1  3.094e-10 ***
## dispersal      0.0146  1      0.904    
## day:infusion  58.5469  4  5.858e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: simp.norm
#                 Chisq Df Pr(>Chisq)    
# day           62.7925  2  2.316e-14 ***
# infusion     274.8962  2  < 2.2e-16 ***
# temperature   47.2450  1  6.264e-12 ***
# dispersal      0.0397  1      0.842    
# day:infusion  53.2218  4  7.658e-11 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

wsimp1.int.notem <- glmmTMB(simp.norm~day*infusion+dispersal+(1|mesocosm),data=df.div.mw)
anova(wsimp1.int,wsimp1.int.notem) #sig***
## Data: df.div.mw
## Models:
## wsimp1.int.notem: simp.norm ~ day * infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1.int.notem 12 378.91 419.13 -177.46   354.91                             
## wsimp1.int       13 349.46 393.04 -161.73   323.46 31.451      1  2.046e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wsimp1.int.nodis <- glmmTMB(simp.norm~day*infusion+temperature+(1|mesocosm),data=df.div.mw)
anova(wsimp1.int,wsimp1.int.nodis) #0.84
## Data: df.div.mw
## Models:
## wsimp1.int.nodis: simp.norm ~ day * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1.int.nodis 12 347.48 387.70 -161.74   323.48                         
## wsimp1.int       13 349.46 393.04 -161.73   323.46 0.0146      1      0.904
##all of the below was worse before normalizing
shapiro.test(residuals(wsimp1.int))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wsimp1.int)
## W = 0.93884, p-value = 9.519e-08
qqnorm(residuals(wsimp1.int))
qqline(residuals(wsimp1.int), col="red")

wsimp1.int.resid <- simulateResiduals(fittedModel = wsimp1.int, plot = T) 

##posthoc things
wsimp1.em <- emmeans(wsimp1.int,~infusion+temperature)
multcomp::cld(wsimp1.em)
##  infusion temperature  emmean     SE  df lower.CL upper.CL .group
##  SG       C           -1.2595 0.0939 198   -1.445  -1.0743  1    
##  SG       H           -0.6695 0.0943 198   -0.855  -0.4835   2   
##  OL       C           -0.0999 0.0934 198   -0.284   0.0842    3  
##  SW       C            0.4404 0.0934 198    0.256   0.6245     4 
##  OL       H            0.4901 0.0937 198    0.305   0.6750     4 
##  SW       H            1.0304 0.0937 198    0.846   1.2153      5
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
 # infusion temperature emmean     SE  df lower.CL upper.CL .group
 # SG       C           -1.300 0.0886 198   -1.474  -1.1249  1    
 # SG       H           -0.692 0.0890 198   -0.867  -0.5162   2   
 # OL       C           -0.105 0.0881 198   -0.278   0.0691    3  
 # SW       C            0.463 0.0881 198    0.289   0.6365     4 
 # OL       H            0.503 0.0884 198    0.329   0.6777     4 
 # SW       H            1.071 0.0884 198    0.896   1.2452      5

3.3.2.2 Water Simpson’s stats - rarefied

str(df.rare.div.mw) #make sure day is a factor still
## 'data.frame':    211 obs. of  27 variables:
##  $ sample_name      : chr  "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
##  $ Observed         : num  25 27 22 28 31 31 28 29 27 21 ...
##  $ Shannon          : num  1.64 1.51 1.56 1.5 1.77 ...
##  $ InvSimpson       : num  2.69 2.36 2.41 2.27 2.98 ...
##  $ orgname          : chr  "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
##  $ mesocosm         : chr  "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
##  $ type             : chr  "Microbial Water" "Microbial Water" "Microbial Water" "Microbial Water" ...
##  $ stage            : chr  "Water" "Water" "Water" "Water" ...
##  $ sex              : chr  "" "" "" "" ...
##  $ date.collected   : chr  "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
##  $ day              : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
##  $ abdomen.length.mm: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ wing.length      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ hindleg.length   : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ special          : chr  "" "" "" "" ...
##  $ dispersal        : chr  "N" "N" "D" "D" ...
##  $ temperature      : chr  "C" "C" "C" "C" ...
##  $ infusion         : chr  "OL" "OL" "OL" "OL" ...
##  $ raw_miseq_reads  : int  30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
##  $ X16scq           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct1         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Wolb_ct2         : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Act_ct           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ notes            : chr  "" "" "" "" ...
##  $ newday           : chr  "early" "early" "early" "early" ...
##  $ lib_size         : num  25797 34139 15861 49298 44559 ...
##  $ is.neg           : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##normalizing because residuals bad below
bestNormalize(df.rare.div.mw$InvSimpson) #orderNorm
## Best Normalizing transformation with 211 Observations
##  Estimated Normality Statistics (Pearson P / df, lower => more normal):
##  - arcsinh(x): 1.697
##  - Box-Cox: 1.6111
##  - Center+scale: 1.5965
##  - Double Reversed Log_b(x+a): 1.5648
##  - Exp(x): 3.1941
##  - Log_b(x+a): 1.6639
##  - orderNorm (ORQ): 1.3294
##  - sqrt(x + a): 1.5573
##  - Yeo-Johnson: 1.5903
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##  
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties 
##  - Original quantiles:
##    0%   25%   50%   75%  100% 
## 1.367 1.957 2.684 3.244 5.127
df.rare.div.mw$simp.norm <- bestNormalize(df.rare.div.mw$InvSimpson)$x.t
shapiro.test(df.rare.div.mw$simp.norm)
## 
##  Shapiro-Wilk normality test
## 
## data:  df.rare.div.mw$simp.norm
## W = 0.99978, p-value = 1
#hist(df.rare.div.mw$Observed)
#shapiro.test(log(df.rare.div.mw$Observed))
#shapiro.test(df.rare.div.mw$Observed)

##full model
wsimp1<-glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
##no dispersal
wsimp1nd<-glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm), data=df.rare.div.mw)
##minus infusion
wsimp1ni<-glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus temperature 
wsimp1nt<-glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)
##minus time 
wsimp1nnd<-glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mw)

anova(wsimp1,wsimp1nd) #0.87
## Data: df.rare.div.mw
## Models:
## wsimp1nd: simp.norm ~ day + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1nd  8 385.60 412.41 -184.80   369.60                         
## wsimp1    9 387.57 417.74 -184.79   369.57 0.0268      1     0.8699
anova(wsimp1,wsimp1ni) #2.2e-16***
## Data: df.rare.div.mw
## Models:
## wsimp1ni: simp.norm ~ day + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1ni  7 487.79 511.25 -236.89   473.79                             
## wsimp1    9 387.57 417.74 -184.79   369.57 104.21      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wsimp1,wsimp1nt) #8.56e-09***
## Data: df.rare.div.mw
## Models:
## wsimp1nt: simp.norm ~ day + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1nt  8 418.72 445.53 -201.36   402.72                             
## wsimp1    9 387.57 417.74 -184.79   369.57 33.144      1   8.56e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wsimp1,wsimp1nnd) #3.65e-10***
## Data: df.rare.div.mw
## Models:
## wsimp1nnd: simp.norm ~ temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1nnd  7 427.04 450.50 -206.52   413.04                             
## wsimp1     9 387.57 417.74 -184.79   369.57 43.464      2  3.646e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(wsimp1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                Chisq Df Pr(>Chisq)    
## day          50.7854  2  9.377e-12 ***
## temperature  42.1861  1  8.299e-11 ***
## infusion    231.0011  2  < 2.2e-16 ***
## dispersal     0.0268  1       0.87    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#                Chisq Df Pr(>Chisq)    
# day          50.7856  2  9.377e-12 ***
# temperature  42.1866  1  8.297e-11 ***
# infusion    231.0025  2  < 2.2e-16 ***
# dispersal     0.0268  1       0.87    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##interactions?
wsimp1.int <- glmmTMB(simp.norm~day*infusion+temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int1 <- glmmTMB(simp.norm~day+infusion*temperature+dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day+infusion+temperature*dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int3 <- glmmTMB(simp.norm~day*temperature+infusion+dispersal+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int4 <- glmmTMB(simp.norm~day*dispersal+temperature+infusion+(1|mesocosm),data=df.rare.div.mw)
wsimp1.int5 <- glmmTMB(simp.norm~day+dispersal*infusion+temperature+(1|mesocosm),data=df.rare.div.mw)

anova(wsimp1,wsimp1.int) #1.148e-09***
## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int: simp.norm ~ day * infusion + temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## wsimp1      9 387.57 417.74 -184.79   369.57                             
## wsimp1.int 13 347.98 391.56 -160.99   321.98 47.592      4  1.148e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(wsimp1,wsimp1.int1) #ns
## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int1: simp.norm ~ day + infusion * temperature + dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 387.57 417.74 -184.79   369.57                         
## wsimp1.int1 11 387.64 424.51 -182.82   365.64 3.9376      2     0.1396
anova(wsimp1,wsimp1.int2) #ns
## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int2: simp.norm ~ day + infusion + temperature * dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 387.57 417.74 -184.79   369.57                         
## wsimp1.int2 10 388.86 422.38 -184.43   368.86 0.7141      1     0.3981
anova(wsimp1,wsimp1.int3) #ns
## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int3: simp.norm ~ day * temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 387.57 417.74 -184.79   369.57                         
## wsimp1.int3 11 391.02 427.89 -184.51   369.02 0.5552      2     0.7576
anova(wsimp1,wsimp1.int4) #ns
## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int4: simp.norm ~ day * dispersal + temperature + infusion + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 387.57 417.74 -184.79   369.57                         
## wsimp1.int4 11 389.79 426.66 -183.90   367.79 1.7812      2     0.4104
anova(wsimp1,wsimp1.int5) #ns
## Data: df.rare.div.mw
## Models:
## wsimp1: simp.norm ~ day + temperature + infusion + dispersal + (1 | mesocosm), zi=~0, disp=~1
## wsimp1.int5: simp.norm ~ day + dispersal * infusion + temperature + (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## wsimp1       9 387.57 417.74 -184.79   369.57                         
## wsimp1.int5 11 391.00 427.87 -184.50   369.00 0.5773      2     0.7493
Anova(wsimp1.int)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##                 Chisq Df Pr(>Chisq)    
## day           71.6034  2  2.828e-16 ***
## infusion     229.6298  2  < 2.2e-16 ***
## temperature   40.5261  1  1.940e-10 ***
## dispersal      0.0345  1     0.8526    
## day:infusion  57.0587  4  1.203e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: simp.norm
#                 Chisq Df Pr(>Chisq)    
# day           71.6032  2  2.829e-16 ***
# infusion     229.6270  2  < 2.2e-16 ***
# temperature   40.5264  1  1.940e-10 ***
# dispersal      0.0345  1     0.8526    
# day:infusion  57.0591  4  1.202e-11 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#tapply(df.rare.div.mw$log(simp.norm), df.rare.div.mw$infusion, mean)
#tapply(df.rare.div.mw$log(simp.norm), df.rare.div.mw$day, mean)

##model checking
shapiro.test(residuals(wsimp1.int)) #awesome
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(wsimp1.int)
## W = 0.93335, p-value = 3.21e-08
qqnorm(residuals(wsimp1.int))
qqline(residuals(wsimp1.int), col="red")

plotResiduals(wsimp1.int)

wsimp1.int.resid <- simulateResiduals(fittedModel = wsimp1.int, plot = T)

##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wsimp1.int.em <- emmeans(wsimp1.int,~infusion+temperature)

multcomp::cld(wsimp1.int.em)
##  infusion temperature  emmean     SE  df lower.CL upper.CL .group
##  SG       C           -1.2646 0.0937 198   -1.449  -1.0798  1    
##  SG       H           -0.6694 0.0941 198   -0.855  -0.4839   2   
##  OL       C           -0.0993 0.0931 198   -0.283   0.0843    3  
##  SW       C            0.4376 0.0931 198    0.254   0.6213     4 
##  OL       H            0.4959 0.0935 198    0.312   0.6803     4 
##  SW       H            1.0328 0.0935 198    0.848   1.2172      5
## 
## Results are averaged over the levels of: day, dispersal 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.

3.3.2.3 Mosquito Simpson’s

From manuscript pre-revisions: “In contrast, the Simpson’s index, which integrates evenness into a measure of alpha diversity, of the adult mosquito microbiome was lower in male mosquitoes relative to female mosquitoes (p<0.0001). This trend was associated with the increased dominance of wAlbB in males. The Simpson’s index of the adult mosquito microbiome did not vary significantly with the dispersal, aquatic chemistry, temperature, or time of emergence.”

hist(df.div.mq$InvSimpson)

df.div.mq$simp.norm <- bestNormalize(df.div.mq$InvSimpson)$x.t #orderNorm
hist(df.div.mq$simp.norm)

##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+(1|mesocosm), data=df.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)

anova(msimp.all,msimp.nodis) #0.10
## Data: df.div.mq
## Models:
## msimp.nodis: simp.norm ~ newday + temperature + infusion + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## msimp.nodis  9 507.26 536.72 -244.63   489.26                         
## msimp.all   10 506.57 539.30 -243.28   486.57 2.6991      1     0.1004
anova(msimp.all,msimp.noinf) #0.26
## Data: df.div.mq
## Models:
## msimp.noinf: simp.norm ~ newday + temperature + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## msimp.noinf  8 505.29 531.47 -244.65   489.29                         
## msimp.all   10 506.57 539.30 -243.28   486.57 2.7247      2     0.2561
anova(msimp.all,msimp.nosex) #2.72e-14***
## Data: df.div.mq
## Models:
## msimp.nosex: simp.norm ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## msimp.nosex:     mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## msimp.nosex  9 562.49 591.95 -272.25   544.49                             
## msimp.all   10 506.57 539.30 -243.28   486.57 57.925      1  2.723e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(msimp.all,msimp.notem) #0.58
## Data: df.div.mq
## Models:
## msimp.notem: simp.norm ~ newday + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.notem  9 504.87 534.33 -243.44   486.87                        
## msimp.all   10 506.57 539.30 -243.28   486.57 0.306      1     0.5801
anova(msimp.all,msimp.notim) #0.31
## Data: df.div.mq
## Models:
## msimp.notim: simp.norm ~ temperature + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## msimp.notim  8 504.93 531.12 -244.47   488.93                         
## msimp.all   10 506.57 539.30 -243.28   486.57 2.3679      2     0.3061
Anova(msimp.all)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##               Chisq Df Pr(>Chisq)    
## newday       2.3823  2    0.30387    
## temperature  0.3062  1    0.58000    
## infusion     2.7438  2    0.25363    
## dispersal    2.8697  1    0.09026 .  
## sex         67.4478  1    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Response: simp.norm
#               Chisq Df Pr(>Chisq)    
# newday       2.3823  2    0.30387    
# temperature  0.3062  1    0.58000    
# infusion     2.7438  2    0.25363    
# dispersal    2.8697  1    0.09026 .  
# sex         67.4477  1    < 2e-16 ***

#AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)
#no time & no temp tied for best, then no infusion

##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(msimp.all)
## W = 0.98791, p-value = 0.0961
qqnorm(residuals(msimp.all))
qqline(residuals(msimp.all), col="red")

msimp.all.resid <- simulateResiduals(fittedModel = msimp.all, plot = T) 

##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)
##  infusion temperature  emmean    SE  df lower.CL upper.CL .group
##  SG       H           -0.2061 0.152 185   -0.506   0.0939  1    
##  SG       C           -0.1150 0.167 185   -0.445   0.2150  1    
##  OL       H           -0.0782 0.143 185   -0.361   0.2045  1    
##  OL       C            0.0129 0.115 185   -0.215   0.2406  1    
##  SW       H            0.1484 0.162 185   -0.172   0.4682  1    
##  SW       C            0.2394 0.193 185   -0.141   0.6197  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
#no diffs

3.3.2.4 Mosquito Simpson’s rarefied

hist(df.rare.div.mq$InvSimpson)

df.rare.div.mq$simp.norm <- bestNormalize(df.rare.div.mq$InvSimpson)$x.t #orderNorm
hist(df.rare.div.mq$simp.norm)

##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+(1|mesocosm), data=df.rare.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.rare.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.rare.div.mq)

anova(msimp.all,msimp.nodis) #0.095 .
## Data: df.rare.div.mq
## Models:
## msimp.nodis: simp.norm ~ newday + temperature + infusion + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## msimp.nodis  9 506.90 536.36 -244.45   488.90                           
## msimp.all   10 506.12 538.85 -243.06   486.12 2.7861      1    0.09508 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(msimp.all,msimp.noinf) #0.25
## Data: df.rare.div.mq
## Models:
## msimp.noinf: simp.norm ~ newday + temperature + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## msimp.noinf  8 504.89 531.07 -244.44   488.89                         
## msimp.all   10 506.12 538.85 -243.06   486.12 2.7719      2     0.2501
anova(msimp.all,msimp.nosex) #2.21e-14***
## Data: df.rare.div.mq
## Models:
## msimp.nosex: simp.norm ~ newday + temperature + infusion + dispersal + (1 | , zi=~0, disp=~1
## msimp.nosex:     mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## msimp.nosex  9 562.45 591.91 -272.23   544.45                             
## msimp.all   10 506.12 538.85 -243.06   486.12 58.334      1  2.212e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(msimp.all,msimp.notem) #0.56
## Data: df.rare.div.mq
## Models:
## msimp.notem: simp.norm ~ newday + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## msimp.notem  9 504.45 533.91 -243.22   486.45                        
## msimp.all   10 506.12 538.85 -243.06   486.12 0.333      1     0.5639
anova(msimp.all,msimp.notim) #0.31
## Data: df.rare.div.mq
## Models:
## msimp.notim: simp.norm ~ temperature + infusion + dispersal + sex + (1 | mesocosm), zi=~0, disp=~1
## msimp.all: simp.norm ~ newday + temperature + infusion + dispersal + sex + , zi=~0, disp=~1
## msimp.all:     (1 | mesocosm), zi=~0, disp=~1
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## msimp.notim  8 504.45 530.64 -244.23   488.45                         
## msimp.all   10 506.12 538.85 -243.06   486.12 2.3387      2     0.3106
Anova(msimp.all)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: simp.norm
##               Chisq Df Pr(>Chisq)    
## newday       2.3527  2    0.30840    
## temperature  0.3332  1    0.56376    
## infusion     2.7917  2    0.24762    
## dispersal    2.9784  1    0.08438 .  
## sex         67.9982  1    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#               Chisq Df Pr(>Chisq)    
# newday       2.3527  2    0.30840    
# temperature  0.3332  1    0.56376    
# infusion     2.7917  2    0.24762    
# dispersal    2.9784  1    0.08438 .  
# sex         67.9982  1    < 2e-16 ***

#AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)
#no time & no temp tied for best, then no infusion

##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(msimp.all)
## W = 0.98823, p-value = 0.107
qqnorm(residuals(msimp.all))
qqline(residuals(msimp.all), col="red")

msimp.all.resid <- simulateResiduals(fittedModel = msimp.all, plot = T) 

##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)
##  infusion temperature  emmean    SE  df lower.CL upper.CL .group
##  SG       H           -0.2117 0.152 185   -0.511   0.0879  1    
##  SG       C           -0.1168 0.167 185   -0.446   0.2128  1    
##  OL       H           -0.0779 0.143 185   -0.360   0.2045  1    
##  OL       C            0.0170 0.115 185   -0.210   0.2445  1    
##  SW       H            0.1446 0.162 185   -0.175   0.4641  1    
##  SW       C            0.2395 0.193 185   -0.140   0.6193  1    
## 
## Results are averaged over the levels of: newday, dispersal, sex 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## significance level used: alpha = 0.05 
## NOTE: If two or more means share the same grouping symbol,
##       then we cannot show them to be different.
##       But we also did not show them to be the same.
#no diffs

3.4 Infusion water things

df.div.inf <- subset(df.div,type=="Infusion water")

df.div.inf.se <- summarySE(df.div.inf,measurevar="Observed",groupvars=c("infusion"))
df.div.inf.se
##   infusion N  Observed       sd        se        ci
## 1       OL 3 45.666667 3.214550 1.8559215  7.985386
## 2       SG 3 23.666667 6.350853 3.6666667 15.776393
## 3       SW 3  4.333333 1.154701 0.6666667  2.868435
# gg.iw.rich <- ggplot(df.div.inf,aes(x=infusion,y=Observed,fill=infusion))+
#   geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#   geom_errorbar(data=df.div.inf.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#   geom_point(data=df.div.inf.se,size=2.5,position=position_dodge(width=0.6))+
#   #geom_boxplot(outlier.shape=NA,alpha=0.5)+
#   theme_cowplot()+
#   ggtitle("Infusion water")
# gg.iw.rich

4 Diversity correlations

4.1 Diversity x day

Doesn’t look interesting

ggplot(df.div.mq,aes(x=as.numeric(day),y=Observed))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(df.div.mq,aes(x=as.numeric(day),y=InvSimpson))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

4.2 Diversity - mesocosm vs. mosquito

library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:cowplot':
## 
##     stamp
## The following objects are masked from 'package:plyr':
## 
##     rename, round_any
meso.corr <- merge(df.div.mq,df.div.mw,by="mesocosm")
ggplot(meso.corr,aes(x=InvSimpson.x,y=InvSimpson.y))+
  geom_point()+
  #facet_wrap(~day.y)+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(meso.corr,aes(x=Observed.x,y=Observed.y))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

5 Session info

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Pacific/Honolulu
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] reshape_0.8.9       bestNormalize_1.9.1 DHARMa_0.4.6       
##  [4] bbmle_1.0.25.1      emmeans_1.10.3      glmmTMB_1.1.9      
##  [7] phyloseq_1.48.0     cowplot_1.1.3       ggpubr_0.6.0       
## [10] ggplot2_3.5.1       Rmisc_1.5.1         plyr_1.8.9         
## [13] lattice_0.22-6      car_3.1-2           carData_3.0-5      
## 
## loaded via a namespace (and not attached):
##   [1] jsonlite_1.8.8          magrittr_2.0.3          TH.data_1.1-2          
##   [4] estimability_1.5.1      farver_2.1.2            nloptr_2.1.1           
##   [7] rmarkdown_2.27          zlibbioc_1.50.0         vctrs_0.6.5            
##  [10] multtest_2.60.0         minqa_1.2.7             rstatix_0.7.2          
##  [13] butcher_0.3.4           htmltools_0.5.8.1       broom_1.0.6            
##  [16] Rhdf5lib_1.26.0         rhdf5_2.48.0            sass_0.4.9             
##  [19] parallelly_1.38.0       bslib_0.8.0             sandwich_3.1-0         
##  [22] lubridate_1.9.3         zoo_1.8-12              cachem_1.1.0           
##  [25] TMB_1.9.14              igraph_2.0.3            mime_0.12              
##  [28] lifecycle_1.0.4         iterators_1.0.14        pkgconfig_2.0.3        
##  [31] gap_1.5-3               Matrix_1.7-0            R6_2.5.1               
##  [34] fastmap_1.2.0           rbibutils_2.2.16        shiny_1.9.1            
##  [37] GenomeInfoDbData_1.2.12 future_1.34.0           digest_0.6.36          
##  [40] numDeriv_2016.8-1.1     colorspace_2.1-1        S4Vectors_0.42.1       
##  [43] qgam_1.3.4              vegan_2.6-6.1           labeling_0.4.3         
##  [46] timechange_0.3.0        fansi_1.0.6             httr_1.4.7             
##  [49] abind_1.4-5             mgcv_1.9-1              compiler_4.4.1         
##  [52] rngtools_1.5.2          withr_3.0.1             doParallel_1.0.17      
##  [55] backports_1.5.0         highr_0.11              ggsignif_0.6.4         
##  [58] MASS_7.3-61             lava_1.8.0              biomformat_1.32.0      
##  [61] permute_0.9-7           tools_4.4.1             ape_5.8                
##  [64] httpuv_1.6.15           future.apply_1.11.2     nnet_7.3-19            
##  [67] glue_1.7.0              promises_1.3.0          nlme_3.1-165           
##  [70] rhdf5filters_1.16.0     grid_4.4.1              cluster_2.1.6          
##  [73] reshape2_1.4.4          ade4_1.7-22             generics_0.1.3         
##  [76] recipes_1.1.0           gtable_0.3.5            nortest_1.0-4          
##  [79] class_7.3-22            tidyr_1.3.1             data.table_1.15.4      
##  [82] utf8_1.2.4              XVector_0.44.0          BiocGenerics_0.50.0    
##  [85] foreach_1.5.2           pillar_1.9.0            stringr_1.5.1          
##  [88] later_1.3.2             splines_4.4.1           dplyr_1.1.4            
##  [91] survival_3.7-0          tidyselect_1.2.1        Biostrings_2.72.1      
##  [94] knitr_1.48              gridExtra_2.3           IRanges_2.38.1         
##  [97] xfun_0.46               Biobase_2.64.0          hardhat_1.4.0          
## [100] timeDate_4032.109       stringi_1.8.4           UCSC.utils_1.0.0       
## [103] yaml_2.3.10             boot_1.3-30             evaluate_0.24.0        
## [106] codetools_0.2-20        tibble_3.2.1            multcompView_0.1-10    
## [109] cli_3.6.3               rpart_4.1.23            Rdpack_2.6             
## [112] xtable_1.8-4            munsell_0.5.1           jquerylib_0.1.4        
## [115] Rcpp_1.0.13             GenomeInfoDb_1.40.1     globals_0.16.3         
## [118] bdsmatrix_1.3-7         parallel_4.4.1          gower_1.0.1            
## [121] gap.datasets_0.0.6      doRNG_1.8.6             lme4_1.1-35.5          
## [124] listenv_0.9.1           mvtnorm_1.2-5           ipred_0.9-15           
## [127] scales_1.3.0            prodlim_2024.06.25      purrr_1.0.2            
## [130] crayon_1.5.3            rlang_1.1.4             multcomp_1.4-26